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ABSTRACT 

The Rossby wave instability (RWI) is a promising mechanism for producing large-scale 
vortices in protoplanetary discs. The instability operates around a density bump in 
the disc, and the resulting vortices may facilitate planetesimal formation and angular 
momentum transfer in the disc dead zone. Most previous works on the RWI deal 
with two-dimensional (height-integrated) discs. However, vortices may have different 
dynamical behaviours in 3D than in 2D. Recent numerical simulations of the RWI 
in 3D global discs by Meheut et al. have revealed intriguing vertical structure of the 
vortices, including appreciable vertical velocities. In this paper we present a linear 
analysis of the RWI in 3D global models of isothermal discs. We calculate the growth 
rates of the Rossby modes (of various azimuthal wave numbers tti = 2 — 6) trapped 
around the fiducial density bump and carry out 3D numerical simulations to compare 
with our linear results. We show that the 3D RWI growth rates are only slightly 
smaller than the 2D growth rates, and the velocity structures seen in the numerical 
simulations during the linear phase are in agreement with the velocity eigenfunctions 
obtained in our linear calculations. This numerical benchmark shows that numerical 
simulations can accurately describe the instability. The angular moixLcntuixL transfer 
rate associated with Rossby vortices is also studied. 
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1 INTRODUCTION 

Vortices may play an importance role in protoplanetary 
discs and planet formation. First, long-lived anticyclonic 
vortices concentrate dust grains in thei r centers and accel- 
erate the growth o f mctcr- sized solids ([Ba rge & Sommerial 
' l995| : lTanga et aLlllOTG: Bracco et al.lll999l: Godon fc Livio 



A promising mechanism for produci ng vortices is 
the Rossby Wave Insta bility (RWI) (jLovelace et al.l 
Il999l iLi et al.1 I2OO0I . I2OOII ). The RWI is fundamentally 
related to Kelvin-Helmholtz instability of shearing 
flows and Rayleigh 's infle cti on point theorem (e.g. , 
Papaloizou fc Pringld Il985l : iPapaloizou fc LinI Il989l : 



20od : Ijohansen et al.ll2004l : iHeng fc KenvonI [20101 '). leading 



to the formation of planetesimals. Second, vortices may gen- 
erate disordered (turbulent) flow, which induces angular mo- 
mentum transfer in the radial direction. This is pa rticularly 
relevant to accretion in the disc dead zone (e.g. iGammid 
119961 : lTerauemll2008l ) , where the gas is not ionised by stellar 
radiation nor by cosmic rays, and turbulence driven by mag- 
netorotational instability is ineffective. Finally, vortices can 
also form at t he edge of planetary gaps (do Va .l-Borro et al.l 
I2OO7I : lYuet al, 2010 : Lin fc Papaloizou 2011a, . b|). thereby af- 
fecting the rate of planet migration. In these different con- 
texts, but mainly for the concentration of solids particles, it 
is important to understand the 3D velocity structure of the 
vortices. 
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Lithwickll2009f ). For two-dimensional (vertically integrated) 
barotropic discs, the RWI relies on the existence of an 
extremum in the background fluid vortensity, deflned by 

(V X v) ■ 2 _ K^ 



c = 



2nE' 



(1) 



where E is the disc surface density, v the fluid velocity, D, 
the disc rotation rate and 



'^'-'-^iir'^) 



dr 



(2) 



is the square of the radial epicyclic frequency (so that k^ /2Q, 
is the vorticity). Since Rossby waves propagate along the 
gradient of vortensity, the instability can be understood 
as arising from the interaction between two Rossby waves 
standing on each side of the vortensity extremum. In pro- 
toplanetray discs, the vortensity extremum is expected to 
form at the boundaries of the dead zone, as the differ- 
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ence in accretio n rate in these region s forms a bump in the 
density profile jVarniere et alj l2006l : llnaba fc Baxg3 l2006l : 
iTerauemI [200i: Kretke et alj 120091 ^ Recent MHD simula- 
tions ( Dzvurkevich et al.ll2010r ) are beginning to address var- 
ious complex physics related to the formation of this den- 
sity/pressure bump. 

So far most of the studies of the RWI have con- 
sidered t wo-dimensi onal (height-integrated) discs (e.g. 
Lvra et al. 2009: Ta gger fc Melial 12009 : IVarniere fc Taggeij 
20061 : iRegalv et all 1201 ll). 



However, while 2D vortices are 



long-lived (but see IChang fc Oishi| 20ld). the situation is 
not clear in three - dimensions (seelBarranco fc MarcuduOOa : 
IShen et al.l l2006l : iLesur fc Papaloizoul |2009|) , and it is im- 
portant to study the formation of vortic es and their sta - 
bility in 3D discs. On the analytical side, lUmurhanl ()201G| ) 
used local shallow-water approximation to examine the RWI 
in idealized situat ions where the density bump is replaced 
by step functions. iLithwickl ([2009!) considered vertically un- 
stratified discs in shearing box approximation. Some 3D as- 
pects of other instabilities controlled by corotation, such as 
the Papaloizou-Pringle instability in r otating tori, have also 
been studied in previous w o rks (e.g. [ Papaloizou fc Pringld 
1 19851 : iGoldreich erahl Il986l : iLatter fc BalbusI l2009l ). How- 
ever, the 3D aspect of the RWI, where Rossby waves are 
self trapped around the density bump, has not been stud- 
ied. 

Recently, Meheut et al. (2010, 2012) carried out global 
3D numerical simulations of the RWI in a full 3D disc. The 
simulations showed that the RWI can develop in 3D discs 
as in 2D, but the resulting Rossby vortices have significant 
vertical structure, contrary to what was expected. In partic- 
ular, the vertical velocity plays an important role in the long- 
term evolution of the vo rtices (Meheut et al., 2012). Indeed, 
the e lliptical instability ()Kerswel]||2002l : iLesur fc Papaloizoul 
[2OO9I) that is responsible for the decay of unstratified (or 
with no vertical flow) vortices is a 3D mechanism, and may 
be affected by the vertical structure. Moreover, the vertical 
velocity can modify the concentration of solids inside the 
vortices, thereby influencing the formation of planetesimals. 

The goal of this paper is to provide an understanding of 
the 3D structure of Rossby vortices that has been observed 
in the 3D simulations of the RWI. To this end, we carry out 
global linear analysis of the RWI in full 3D discs and com- 
pare our results with those obtained from numerical simu- 
lations. Our paper is organized as follows. We first present 
the 3D linear perturbation equations (Section[2)| and the nu- 
merical methods for solving the eigenvalue problem (Section 
[3]). The results of our linear calculations are described and 
discussed in Section |4] In order to compare with the results 
of Meheut et al. (2010), we present in Section [5] full numeri- 
cal simulations of the development of Rossby vortices using 
the same setup and parameters as our linear analysis. We 
conclude in Section 6. 



2 EQUATIONS 

The linear perturbation equ ations for stratified 3D discs 
are similar to those given in IZhang fc Lail pOOd ) in their 
study of the tidal excitation of 3D wa ves in discs (see also 
lOkazaki et al.lll987l : lTanaka et al.l2002h . Here we summarize 
the notations and key equations relevant to our analysis. 



2.1 Governing equations 

We consider a geometrically thin gas disc and adopt cylin- 
drical coordinates {r,tp,z). The governing equations read: 



O'v 

p-^ + (v • V)pv = -Vp - pV^G 

|+V.(pv)=0, 



(3) 
(4) 



where p is the density, p the pressure, v the velocity. The 
disc is assumed to be non-self-gravitating, the gravitational 
potential $0 is the due to the central star only. We will 
consider an isothermal equation of state, p = c^p, with c^ 
the constant sound speed. 

The unperturbed disc has velocity Ve — (0,rQ,0), 
where the angular velocity Q, — 0,{r) is taken to be a func- 
tion of r alone. 

Thus the vertical density profile is given by 

Pe{r,z)= exp(-Z^/2), with Z = z/h (5) 

\'2-Kh 

where h = h{r) — Cs/0,± is the disc scale height, E = E(r) = 
J dz pe is the surface density, and i^± is the vertical oscil- 
lation frequency of the disc and is equal to the Keplerian 
frequency fix ~ {GM/r'^Y''^ where G is the gravitational 
constant and M the central star mass. 



2.2 Linear pertubations equations 

We now consider linear perturbation of the disc. For sim- 
plicity, we shall assume that the perturbation is isothermal, 
so that the (Eulerian) density and pressure perturbations 
are related by 5p = (?a5p. The linear perturbation equations 
read 

^" + (ve-V)u-f (u-V)ve = -Vr?, (6) 



dt 
dSp 
dt 



+ \/ ■ {p,U + Ve5p) =0, 



(7) 



where u = (5v is the (Eulerian) velocity perturbation, and 
77 = Sp/pe is the enthalpy perturbation. Without loss of 
generality, the perturbation variables are assumed to depend 
on t and ifi as 



UjTj, Sp (X exp{imifi — iuit) , 



(8) 



where the azimuthal mode number m > is an integer and 
uj is the (complex) wave frequency. Equations Q and (0 
then reduce to 

d 



or 



(9) 
(10) 
(11) 



—iuu^ + -—Ur = r], 

zil r 

d 

az 
.- /9e 19, , im d , , „ /i„N 

-ItJ— -77-1- -—-{rpEUr) H peU^ + —-[PeU^) = 0. (12) 

ci r or r az 

Here uj is the "Doppler-shifted" frequency 

Lu = UJ ~ mQ, (13) 

and K is the radial epicyclic frequency defined by equa- 
tion (2). 

We consider radiative boundary condition such that 
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Figure 1. Radial profile of the surface density and vortensity f . 
The radius is in units of ro , and the surface density is normalised 
to Sq. 



waves propagate away from the density bu mp at both th e 
inner and outer boundaries of the disk (e.g. jYu &: Lill2009r ). 
They are defined as follows. We take 

dUrO,2 ., dSho.2 ,, c, /■,A\ 
= lkUrO.2 , r— ^ IKO/Io,2 , (14) 

dr dr 

we then substitute tliese equations into the disk perturba- 
tion equations, we will arrive at a complex polynomial at the 
boundary that determines the value of k. Then we choose 
the appropriate root of the polynomial so that the group 
velocity is in the outgoing direction. We also implement the 
null radial velocity boundary condition to check how the 
eigenvalues are affected by different boundary conditions. 

To separate out the z-dependence, we e xpand the per- 
turbations with Hermite polyn omials Hn fsee lOkazaki et al.l 
Il987l : lKatdl200ll : lOgilvidboO^ ') 



Urir,z) 
Uv{r,z) 

Uz(r, z) 



= E 



r?n(r) 

Urn{r) 



HniZ), 



J2^Mr)HUZ), 



(15) 



where H'„ — dHn/dZ, and the Hermite polynomial is defined 
hyHn{Z) = (-l)"e^''/2d"(e-^''/2)/rfZ". Note that Ho = 1, 
Hi = Z and H2 = Z'^ -- 1. 

With the expansion in (jlSp . the fluid equations (fI?))- H12p 
become 



„„ d n/i 

-lUUrn — 2ilU^n = ——-Tin H ??„ 

dr r 

, {n + l)(n + 2)n 






-tU)U,„ = 



h 



(16) 

(17) 
(18) 



. . 77,1 f d , ^ nii\ u d 

-lOJ^r H- -- InrL H Um H llr,n-2 + -rUrn 

d \ar r I r dr 



im n 

r h 



where 



H = d\nh/d\nr. 
© 2011 RAS, MNRAS 000, [iHH] 



(19) 
(20) 



We note here that for the isothermal discs we consider, 
l-i ^ 0. Eliminating u^„ and Uzn from eqs. p6|) - (|19l ). we have 



dri„ 



2mn 



D . 



dr 


ruj 










+ - [nvn + {n+l){n + 2)r]n+2] , 


(21) 


dUrn 


= - 


■dln(rE) mK^ ' 


1 

Urn + T^T 


( m^ n \ 

2 + h2 '?" 
\ r^ h^ J 


dr 


dr ' 2r^0j 



luj H . , 

ci r 



where we have deflned 
D = 



~2 

- id 



K — (oj — mfi) . 



(22) 



(23) 



If we only consider the n = term in the expansion (jlSp . 
then Uz = and eqs. (|2ip - (|22p reduce to the dynamical 
equations for 2D discs. However, since n ^ 0, ea. (|16|) -(|19 [) 
are an infinite set of coupled ordinary differential equations 
and the n — component is coupled to the n — 2, 4, • • • 
components, giving rise to nontrivial vertical structure. In 
our linear calculation below, we will truncate the Hermite 
series at ti = 2 and only include the n — and n = 2 
components. 



3 SETUP AND METHOD OF SOLUTION 

We consider the following Gaussian density bump in the 
disc: 



E/Eo = l + xexp[- 



jr - rof 
20-2 



(24) 



where ro and a are the position and the width of the bump. 
We choose the bump parameters 



X = 0.25, a/ro = 0.05, 



(25) 



and normalize ro = 1. The sound speed is given by 
Cs/{roQo) = 0.1, where f2o is the Keplerian orbital frequency 
at ro . The corresponding surface density and vortensity pro- 
files are shown in Figure [l] 

Equations (|16p -(|19 p are transformed into an eigenvalue 
problem and solved using two different methods and numer- 
ical codes: 

Method ( 1): We adapt th e code (MODINT) originally 
developed by iTagger fc PellatI ([1999|) in their study of the 
accretion-ejection instability in 2D magnetized discs. We 
add the equations associated with the n = 2 components 
in the code, but the method is otherwise similar. The ra- 
dial direction is discretised on a grid logarithmically spaced 
with a resolution of 255 points. The density is defined at the 
centre of the grid cells, whereas the velocities are defined at 
the cell boundaries. Although t he problem is now 3D, it is 
simpler than the one studied in lTagger fc PellatI (| 19991 ). as 
there is no magnetic field. This allows us to solve the set 
of equations on the real integration axis, which give more 
accurate eigenvectors. 

Method (2): We use the relaxation method to solve 
this two-po i nt bo undary eigenvalue problem as detailed in 
IPress et al.l (|l992l ) . The computation cost is relatively low in 
this method and we can afford higher resolution. Typically, 
we use 800 uniform grid points in our calculation. 

Since the Rossby waves are concentrated around the 



4 H. Meheut, C. Yu and D. Lai 




Figure 2. The RWI growth rate (in units of Qq) ^^ ^ function of 
the azinuthal mode number m obtained by linear calculations for 
2D {dashed line) and 3D (solid line) discs. 




density bump, we choose for both codes the inner boundary 
at rin/ro = 0.4 and the outer boundary at rout/ro = 1.6. 
The two codes give the same results when the same null 
radial velocity boundary conditions are adopted. 



Figure 4. Angular momentum flux as defined by eq. (28) for the 
n = (upper panel) and n = 2 (lower panel) components of the 
m = 4 mode. The radius is in units of ro and the flux in units 
of SorgCQ. In the lower panel, the dashed line shows Fq, and the 
insert shows a blowup of F2 around the density bump. 



4 RESULTS OF LINEAR CALCULATIONS 

As noted above, we have implemented two types of boundary 
conditions to solve the linear eigenvalue problem. The radia- 
tive boundary condition is preferred, since it corresponds to 
the spontaneous growth of perturbations around the den- 
sity bump without energy input from the regions outside 
the bump. But whatever the boundary conditions, our cal- 
culations show that for the choice of disc parameters and 
computation domain adopted in this paper, the complex 
eigenvalues oj are almost identical. This can be understood 
since only a small amount of wave energy is leaked out of 
the Rossby zone as outgoing (away from the density bump) 
density waves. 

Figure [2] shows the linear growth rate of the 3D Rossby 
mode trapped around the density bump as a function of 
the azimuthal wavenumber m. The 2D result, obtained by 
including only the n = Q terms in the expansion (|15|l . is 
also shown for comparison. We see that the 3D growth rates 
are only slightly smaller than the 2D ones under the same 
conditions (see Section 4.1 below). In all cases, we find that 
the real part of the mode frequency ujr close to mSlo, with 
Ur/mQ.1) = 0.974 for m = 2 and 0.986 for m = 6. 

Figure [3] depicts the eigenf unctions of the ni = 4 mode. 
For each variable, we show the n = Q and n = 2 components, 
except for the vertical velocity which has a nul n = compo- 
nent. The amplitudes of the eigenfunctions are normalised so 
that the maximum value of the radial velocity perturbation 
|Mro| equals unity (this maximum occurs at r ~ ro). Note 
that the Rossby mode is confined around the corotation ra- 
dius (where Ur = mO.), which is close to the density bump. 
But the mode can leak out as spiral density waves inside the 
inner Lindblad resonance (where Lo — mQ. = —k) and outside 
the outer Lindblad reso nance (where oj — m n = k). [See, e.g., 
iTsang fc Laill2008l and lLai fc Tsang||2009l for discussion on 
the wave zones for Rossby waves and density waves.] The 
n = spiral wave (as a function of r) satisfies the WKB 



amplitude relation for rjo (|Tsang fc Laill2008f ). 



\Vo{r)\ oc 



r2E2 



1/4 



and similar relations for Uro and u,po. The n — 2 component 
is driven by the coupling with the n = component. The 
vertical velocity has a small but non-zero amplitude, and as 
this is a pure n = 2 component, its amplitude increases with 
height whereas the other variables (ur, u^ and 5p) are dom- 
inated in the Rossby wave region by the vertically constant 
component [n = Q). This means that the vertical velocity 
component can not be neglected when the flow pattern is 
of importance, as for instance for the study of the concen- 
tration of dust in the vortices. In other words, the thin disc 
approximation do not capture the full structure of the in- 
stability (see also Section 5). 



4.1 Angular momentum flux 

A useful way to understand the RWI is to examine the an- 
gular momentum flux carried by the wave. The time aver- 
aged transfer rate of the ^-component of angular momen- 
tum across a cylinder of radius r is given by (see, e.g. , 
Lvnden-BeU fc Kalnaid 1 19721 : iGoldreich fc Tremaind Il979l : 
Tanaka et al.ll2002[ ) 



F{r) = (r I dz dippe{r,z)ur(r,if,z,t)u,^{r,ip,z,t)\, 

(27) 
where () stands for the time average. Using 
ur{r,e,z,t) = Re[ur„H„{Z)e'^'^'^-^% u^{r,,fi,z,t) = 
'^s [■"■(= n^"(-^)^]' ^^ ^^'^ that the angular momen- 
tum fl ux associated with the n component of the m mode 
is (see lZhangfc Laill2006l ) 



F„{r) = n\Tvr ERe{urnu'^n)- 



(28) 



(26) 



Figure |4] shows the angular momentum fluxes for the 
n = and 2 components of the m = 4 mode. Around the 
density bump, the flux is dominated by the n — compo- 
nent. Thus it is not surprising that the 3D linear growth rate 
is close to the 2D value (see Figure 2). Indeed, the growth 
of the RWI arises from the positive _Fo around the coro- 
tation (close to the density bump): since the perturbation 
inside (outside) the corotation carries negative (positive) an- 
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(a) Radial velocity 



(b) Aziinutlial velocity 




^•<' '■-< /■■■,■■ 

;; 11/' ! 



(c) Dcnsitiy 



(d) Vertical velocity 



Figure 3. Eigenfunctions of the m = 4 mode: The dotted line shows the real part, the dashed line the imaginary part, and the solid line 
the absolute value. The vertical dot-dashed lines give the positions of the inner and outer Lindblad resonances. The radius is in units of 
ro, see also section |4] 



gular momentum, a net outward angular momentum trans- 
fer across the corotation induces mode growth. 

The finite n — 2 wave component arises from its cou- 

J iling to the n = component. As detailed in IZhang fc Lail 
200a ) (see their section 6.2), n ^ 1 waves propagat- 
ing across the corotation are strongly attenuated, and the 
coro tation acts as a sink for waves with n ^ 1 (see 
also iPapaloizou fc Prin"gi3ll985l : iKatol l2003l : iLi et alllJOO^ : 



iLatter fc Balbusll2009l ). This is consistent with the behavior 
of F2 around the corotation: we see from the insert of Fig- 
ure |4] that F2 changes sign from positive to negative across 
the corotation radius. This explains our numerical result (see 
Figure 2) that the 3D RWI growth rate is smaller than the 
2D growth rate. 



5 COMPARISON TO NUMERICAL 
SIMULATIONS 

We have performed full numerical simulations with the same 
disc configuration and parameters as those used in our linear 
calculations. The comparison to simulations is important for 
different reasons: 



(i) In our linear calculations, we included only the first 
two elements (n = 0, 2) of the Hermite polynomial decom- 
position, the higher order terms being neglected. The non- 
linear simulations can confirm if this approximation is cor- 
rect. 

(ii) The previous simulations of (Meheut et al., 2010) 
showed unexpected vertical structure in the Rossby waves. 
New simulations with the same setup as the linear analysis 
can confirm if the vertical structure is correctly handled and 
if the simulations are not altered by the numerical methods 
or boundary conditions. 



5.1 Numerical methods 

The numerical methods used for this study are very similar 
to the one of Meheut et al., 2012. We used MPl-AMRVAC 
(jKeppens et al.ll2012l ) with the total variation diminishing 
Lax-Friedrich scheme with a third order accurate Koren 
limiter. The grid is cylindrical with r/ro in the range of 
[0.4,1.6], z/ro in the range of [0,0.6], and full azimuthal 
range [0, 27r] is covered. This vertical extension of the grid 
corresponds to 6 scale height at ro. A high vertical extension 
of the grid is needed to resolve the vertical structure of the 
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disc at the outer edge of the grid. The grid is fixed during 
the simulation but the initial resolution is determined to fit 
the regions of high density gradient. The base resolution is 
(128, 32, 32) and up to four level of refinement are allowed, 
reaching a resolution of (2048, 512, 512) near tq. The bound- 
ary conditions are the same as Meheut et al. (2012) with a 
null radial velocity at the boundaries. This boundary con- 
dition is slightly different from the radiative boundary con- 
dition. However, as discussed in Section 4, the two different 
boundary conditions produce nearly identical linear mode 
frequency and growth rate. 

In our simulations, the initial midplane density profile 
is given by 



PMij-) 



2TTh 



1 + X exp 



{r - ro) = 



2cr2 



(29) 



The vertical profile of the density is then chosen to achieve 
hydrostatic equilibrium: 



pe{r,z) = Pa/ (r) exp 



GM / 1 



(30) 



This gives a surface density profile very similar to eq. 24. 
The azimuthal velocity is determined by force balance in 
the radial direction. 

The small difference from the linear analysis resides in 
the upper region of the simulation where the density reaches 
very low values. In the simulations, there is a floor value for 
the density. To check if this difference is of importance we 
have performed some tests, one with a floor density of one 
tenth of its initial value and one where the height of the 
numerical box has been doubled has been doubted. We find 
a negligible difference in the resulting growth rates. 

One simulation has been performed without any pertur- 
bation to check that the initial state is indeed an equilibrium 
state. For the other simulations we add small perturbations 
in the radial velocity: 



ro^o 



esm{mip) exp 



{r - ro)^ 



2cr2 



(31) 



with e ~ 5 X 10"'' and m = 2, 3, • • ■ ,6. 



5.2 Results of the simulations 

Starting from the initial perturbation (31) for a given m, 
the simulation follows the exponential growth of the insta- 
bility and the saturation phase on timescale of a few tens of 
rotations. The exact timescale to reach saturation depends 
on the azimuthal mode number m. A fit of the exponential 
growth can give the growth rate 7 of the RWI. An example 
(for m = 3) is shown in Figure [5] and the density and radial 
velocity perturbation during the linear phase (ffio ~ 18) are 
shown in Figure (6] The azimuthal mode number {m — 3) is 
clearly seen in this figure. 

Figure \7\ shows the linear RWI growth rates obtained 
from our 3D simulations for different m's. We point out 
here that the estimated growth rate depends on the defined 
position of the beginning and the end of the linear phase. 
The growth rates of the different tti's are then given with 
an uncertainty of 5 to 10% due to the difficulty to find the 
best fit for the linear growth. The results are compared with 
the mode growth rates from our linear perturbation calcu- 
lations. We see that the numerical growth rates are in good 



O 




10 20 



30 



40 50 60 



Figure 5. The amplitude of the density perturbation (on a log- 
arithmic scale) as a function of time (in units of l/f2o) in the 
m = 3 simulation. The y-ax\s shows the logarithm of the surface 
density perturbation (5S, where ST, = maXr,(p|S — {T,)^\. The two 
solid lines represent linear fits to the exponential growth phase of 
the RWI. The difEcrencc in the fitted linear growth rates is about 
5%. 




Figure 6. Perturbed density (left) and radial velocity (right) in 
the midplane at tQo ~ 18 in the m = 3 simulation. 



agrement with the linear perturbation theory results. For 
the m ^ 4, the non-linear simulations give slightly smaller 
growth rates as compared to the linear perturbation results. 
This could be due to the higher numerical viscosity on the 
smaller-scale modes or to the higher (n > 2) terms in the 
Hermite decomposition that have not been considered in the 
linear approach that would be more important for higher m. 
Some velocity streamlines in a vertical frame of the disc 
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Figure 9. Numerical 'eigenfunctions' of the m, = 4 mode obtained with the non-linear code at tfio ~ 18: The dotted line is the real part, 
the red dashed line is the imaginary part, and the solid line is the amplitude of the m, = 4 element of the azimuthal Fourier transform 
of the physical quantities. 




Figure 7. Growth rate obtained for different azimuthal modes 
with the 2 methods: the linear analysis {solid line) and the non- 
linear simulations {dashed line). 



are plotted in Figure [S] showing one of the vertical roll struc- 
tures we are studying in this paper. 

To compare the flow velocity and density structures ob- 
tained from our 3D simulations with the linear eigenfunc- 
tions, we can apply Fourier transform in the azimuthal direc- 
tion to our numerical flow outputs. Figure [5] gives an exam- 
ple (for the m = 4 mode at time Qot — 18) of these numer- 
ically determined "eigenfunctions". Note that the density, 



radial and azimuthal velocities are evaluated at the disc mid- 
plane, whereas the vertical velocity is plotted slightly above 
the midplane. This figure should be compare with Figure |3l 
where the same eigenfunctions from the linear perturbation 
theory are shown. Here we have not separated the differ- 
ent vertical components but one can see that the n — 
component dominates except obviously for u^. On Figure 
1101 we have plotted the amplitude of the azimuthal veloc- 
ity eigenfunction obtained with the two methods and with 
the same normalisation, so they can be directly compared. 
However for the linear calculation only n = is taken into 
account whereas the multiple components of the vertical 
structure have not been separated in the numerical simu- 
lations. The good agreement between our numerical simula- 
tions and linear analysis indicates that the simulations can 
correctly describe the RWI. The complicated radial struc- 
ture of the mode in the corotation region coupled with the 
vertical structure is responsible for the complexity of the 
flow, as plotted in Figure [S] 



6 CONCLUSION 

In this paper we have carried out linear analysis of the 
Rossby Wave Instability (RWI) in 3D stratified isothermal 
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Figure 8. Velocity streamlines in a vertical frame at tf7o ~ 18. 



growth rate is due to the absorption of the wave component 
with vertical structure at the corotation resonance. 

We have also calculated the angular momentum flux 
carried by the Rossby vortices. This shows that the growth of 
the instability is due to the exchange of angular momentum 
between two Rossby waves on each side of the density bump. 
The angular momentm transfer tends to reduce the initial 
density bump. This is an important feature of the RWI as 
it may lead to the transfer of angular momentum through 
the dead zone of protoplanetary discs. To fully understand 
this process and the competition with the bump formation 
process, one will need a self-consistent simulation including a 
dead-zone inside a ionised disc. In the mean time, analytical 
model of these processes as the one presented here, can give 
a first insight of this dynamical evolution. 
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discs. Our linear calculations show that the vertical velocity 
inside the Rossby vor tices, that was obtained in previous nu- 
merical simulations (jMeheut et al.ll2010l ) is expected when 
the disc scale height varies with the disc radius. This vertical 
velocity is of crucial importance for understanding the con- 
centration of dust inside the vortices and the formation of 
planetesimals. Detailed comparison with 3D numerical sim- 
ulations show that the simulations of Meheut et al. (2010, 
2012) can correctly handle the RWI in 3D. Hence these sim- 
ulations are robust tools to study the non-linear and long- 
term evolution of this instability. We have also shown that 
the linear growth of Rossby waves is only slightly affected 
by the disc vertical stratification, with similar mode growth 
rates in 3D as in 2D. The small reduction of the 3D mode 
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